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The main aim of tlie paper is to incorporate tlie non-linear kinetic term into non-Markovian 
transport equations described by a continuous time random walk (CTRW) with non-exponential 
waiting time distributions. We consider three different CTRW models with reactions. We derive 
new non-linear Master equations for the mesoscopic density of reacting particles corresponding to 
CTRW with arbitrary jump and waiting time distributions. We apply these equations to the problem 
, of front propagation in the reaction-transport systems with Kolmogorov-Petrovskii-Piskunov (KPP) 

■ kinetics and anomalous diffusion. We have found an explicit expression for the speed of a propagating 

^ ' front in the case of subdiffusive transport. 
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I. INTRODUCTION 

^ ' This article addresses the problem of the mesoscopic description of reaction-transport system of particles performing 
a continuous time random walk (CTRW) [ij. One of the main challenges is an implementation of the description 

^ • of chemical reactions in non-Markovian transport processes governed by CTRW with, non-exponential waiting time 

^ distributions. There exist several approaches and techniques to deal with this problem |2h11|. In particular, there 

-j-J r are many efforts to incorporate the chemical reactions into subdiffusive transport. Different models lead to various 

' fractional reaction-diffusion equations corresponding to kinetic regimes [12 20] and subdifFusion-limited reactions 21 1. 

"j^ , Our main objective here is to discuss how to incorporate the non-linear kinetic term into non-Markovian transport 

' equations which is still an open problem. We consider a one-component reaction-transport system consisting of 

I ' independent particles X that follow CTRW. Let p (x, t) represent the density of these particles at point x and time t. 

'■^ ■ The main purpose is to derive the non-linear Master equation for the density p{x, t) in the following form 
C • 

' where the non-linear evolution operator L has to be determined. The challenge is to derive the Master equation for 

^ . an arbitrary CTRW model coupled with non-linear reaction. We assume that the chemical reaction follows the mass 

' action law and the reaction term is of the form r [p] p. The density of other species participating in the reactions are 

CO . held constant. It has been shown recently that for non-Markovian transport we cannot just add the term r (p) p to 

C ' the right hand side of the evolution equation ^ [l^-O- It is convenient to represent the non- linear reaction rate 

1 r [p) as the difference between the birth rate r+ (p) and the death rate r^ [p) 

a\ ■ 

O : r{p) = r+ (p) - r- {p) . (2) 

a\ . 

• As an illustration, let us consider the classical Schlogl first model [24I 



A + X^A + 2X, x'^B, (3) 

b : 

5^ I where /ci, fc2 and are the reaction rate constants. If we denote the densities of particles X and A hy p and pA 
respectively, then the birth and death rates are 

r+ (p) kipA, (p) = fcs + fc2p. (4) 

In what follows the density pA of the catalyst A is assumed to be constant. 

In what follows we consider three different non-Markovian models for reaction and transport processes. W e apply 
these models to the problem of propagating fronts in reaction-transport systems with non-standard diffusion 1231 . |2J| 
(see also a recent review [1^). The theory of subdiffusive propagation of a front is presented in [1, [ll|, [26l430j |. the 
superdiffusive propagation is studied in [3l|. 
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II. MODEL A 



Non-Markovian behavior of particles performing CTRW occurs when diffusive particles get trapped for random 
times with non-exponential distribution. Let 4> (t) and w (z) denote the waiting time probability density function 
and the dispersal kernel respectively. For simplicity we consider the uncoupled case when jumps and waiting times 
are independent. If the distribution of waiting times is exponential: (pit) = Aexp(— At), the transport model is 
Markovian, and the difficulty of implementation of a non-linear reaction term does not arise. In this case the Master 
equation ([ij takes the form of the Kolmogorov- Feller equation with reaction 

dp f 

— = A / p{x — z, t)w (z) dz — \p + r (p) p (5) 

ot Im 



with a clear separation of the transport and reactive terms on the RHS [23| . In the non-Markovian case these terms are 
not additive (l^. [l5l - [l7| . The key question is how the chemical reaction influences the transport process. For Model 
A, we assume that transport processes associated with CTRW and chemical reactions are independent. This case has 
been considered in the series of papers [l^ [TA .IS] . The main assumption here is that when particles are trapped, the 
waiting time is the same for all particles including newborn particles. One can think of biological applications when 
cells or bacteria are trapped in some confined region, say, at time t < t, they proliferate over the trapping (waiting) 
time t — T, and then they are released at time t. In particular, for the problem of virus infection and its propagation, 
the random waiting time occurs due to virus reproduction inside infected cells [s^ . 



A. Conservation laws for particles 



The first step in the derivation of the non-Markovian Master equation ([T|) is to formulate the integral balance 
equations for the density p{x, t) and the auxiliary density j {x, t) . The latter describes the number of particles arriving 
at point X exactly at time t (see, for example, [ICt Il2l|). The balance equations for p{x,t) and j {x,t) can be written 
in the following form 

p (X, t) = po (x) e^o Hp(-.u))du^ ^ /■* ■ ^) g/J r(p(x,u))du^ [t - t) dr, (6) 

jix,t) = [ pq{x- z)e^o''^Pi^''-))'^-^w{z)(l>{t)dz 
Jr 

+ 11 j{x- z,T)e^>'^p'^'-'-'-'^^'^'^''w{z)(l3(t-T)dzdT. (7) 
JQ Jr 

These equations represent the balance of particles due to non-linear chemical reaction and transport process described 
by CTRW 1:]. Eq. ^ is the conservation law for particles' density p {x, t) at point x at time t. The first term on the 
RHS of ([6]) represents the particles that stay at their initial position up to time t. Their density grows with the rate 
r {p {x,t)) during time interval {0,t). The first term involves also the survival function vj/ (t) — (t) dt which is 
the probability that particles stay at their initial location up to time t. The second term gives the number of particles 
that arrive at point x at previous time t < t and grow with rate r (p) during time interval (r, t) so that no jumps 
take place during this time interval. Note that the initial distribution of particles po (x) is set up in such a way that 
a random walk for all particles starts at i = (no aging effects) . Eq. ([6|) is the conservation law for the particles 
that arrive at point x exactly at time t. The first term on the RHS of ([7]) represents the particles that are at the 
point X — z at time t = 0. Their density increases with the rate r (p) during time interval (0, t) and they jump to the 
point X at time t. The second term describes the particles that arrive at the point x — z at some time t < t and react 
up to time t at which the particles jump to position x. It should be noted that ([6]) and ([7]) are mesoscopic mean-field 
equations. We neglect the internal fluctuations due to the flnite number of particles. In general random fluctuations 
could modify the macroscopic behavior of the reaction-transport systems (see, for example, [s^. Issj). 



B. Non-linear Master equation 



Let us now derive the evolution equation for the density p{x,t). Since the balance equations ([6]) and ^ are non- 
linear, we cannot apply directly the standard technique of Fourier-Laplace transforms. Instead, we differentiate the 
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density p {x, t) given by (|6]) with respect to time 

1^ = J {X, t)+r{p)p^ po (X) e^o r(pM)du^ „ ^' ^ (^^ r(p(.,u))du^ _ ^) (g) 

The last two terms can be interpreted as the density of particles that leave the point x exactly at time t 

i (x, = Po (a;) e-''"' '^(''(^■""''"^ (i) + / J (a;, r) e-''X''(^'"»''"(?!. - r) dr. (9) 

Jo 

It follows from ^ and dH) that 

j {x,t) = / i [x — z,t)w [z) dz. (10) 

Then ([5]) can be rewritten as 

— = / i (x — z,t) w (z) dz — i (x,t) + r (p) p. (11) 
OT Jr 

This equation has a very simple meaning of a balance of particles at point x. The first term on the RHS gives the 
number of particles coming to x from different positions x — z, where the jump size z has the distribution w (z). The 
second term gives the rate at which the particles leave the position x. The last term describes the growth of particles 
due to chemical reactions. In the linear case, a similar equation has been used IT^ [isj . The advantage of having this 
equation is that we do not need the Fourier transform to get the closed equation for the density p{x,t) . We can now 

find an expression for the density i {x, t) in terms of p {x, t). We divide (jG]) and (|9]) by the factor exp ^ r [p (x, u)) duj 

and take the Laplace transform of both equations 



{p (x, t) e- ^0 Kp(-,«))rf«| ^ (j.) ^jr{j {x, t) e' '■(p(^,«))d«|j ^ ^ 



C^^iix,t)e-^o-(p(^.^))dn^^ ^ \^po{x)+c\^j{x,t)e-^o-(pi^.^))'iu^^'^ ^(s), 

where ^ (s) = C{'i> (t)} and (f>{s) — L {0 (i)}. From these equations, we obtain 

0(,s) 



£{z(a;,i)e-^o'-(p(:-,«))*'| 
Inverse Laplace transform gives 



vl'(,s) 



{p(a;,t)e^-''o''(''(^'"»''"|. 



i(x,i)= / X(t-r)p(2;,T)e^>(''("'"»''"dr, 



where K (t) is the standard memory kernel defined by its Laplace transform 

*(s) 1 - <?!>(s) 

Substitution of ^ into (HH) gives us the non-linear Master equation 

^ = / X(<-t) [ /p(x-z,T)e/>(''("-"'"»'^"w(z)d2-p(x,T)e^>(''("'"»'^" 
dt Jo Ur 
+r (p) p. 



(12) 



(13) 



(14) 



This non-linear equation is the main result of this paper. For some particular cases, it can be reduced to known 
equations in the literature. For example, when the reaction rate r (p) — r = const, Eq. (I14p takes the form 



dt In ^ ' 



p{x — z,t)w (z) dz ~ p {x, t) 



dr + rp, 



(15) 
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and this model with constant rate r has been formulated in [l^. Since the effective memory kernel K [t — t) e'"^*"'^^ 
depends on the reaction rate, it is tempting to conclude that this equation describes the coupling of chemical reaction 
and transport. We believe that this conclusion is misleading. In fact this equation describes the complete decoupling 
of transport with memory effects and linear reaction. To show this, let us make a substitution 



p (x, t) = n (x, t) 



Then we obtain the Master equation for n {x, t) 

dn 
'dt 



n{x ~ z,t)w (z) dz — n (x, r) 



dT 



(16) 



(17) 



which is independent of the reaction and describes the transport of passive particles. So we have a complete decoupling 
in which the density p (x, t) is the product of the density of passive particles n (x, t) and the exponential factor e''* 
due to the chemical reaction. 



C. Subdiffusive transport 



Now consider slow anomalous diffusion for which the waiting time pdf (t>{t) has a power-law tail: (j){t) ^ (To/t)^'^'^ 
with 0<7<lasi— >-oo. Clearly, the first moment t(j){t)dt is divergent for < 7 < 1. As a result, the mean square 
displacement < x^ > depends on time t as [H . Here we use the following expression for the survival probability 



^' (t) = E-y 



To 



< 7 < 1, 



where [x] = x^/T {-/n + 1) is the Mittag-LefHer function. The waiting time pdf (j){t) — —-^Ej 
singular behavior V~^a.s t 0. The Laplace transforms of ^E* (t) and (/)(t) are 



To {sTq) 



7-1 



1 + (sToy 



0(s) 



1 



1 + (STO) 



(18) 



has 



(19) 



The advantage of using the Mittag-Leffler function is that we can find the fractional reaction-transport equation 
without passing to the long-time large-distance limit. We find from (jl3p that the Laplace transform of the memory 
kernel is 



K{s) = 



(20) 



Eq. takes the form of a non-linear fractional equation 



dp 
dt 



^ D]-^ ( / p{x-z,t)e-^or(p(^~-'^))dUw{z)dz 



lo r(p{x,u))du 



eJo 



-D 



- lo r{p{x,u))du 



(21) 



where D\ ^ is the Riemann-Liouville fractional derivative defined as 



D\-^p{x,t) 



1 



d 



p{x, T)dT 



r{i~j)dtj, (t-ry 

Now assume that the dispersal kernel w (z) is an even and rapidly decaying function for large z. We expand the 
expression in the brackets (fT4| for small z and truncate the Taylor series at the second moment. We obtain 

(22) 



|J = if (t - r) p (x, r) ef^ '■(p(-^"))^«dr + r (p) p. 



where ct^ — J-^ z^w (z) dz. Note that a similar non-linear equation has been derived in [igI [l8l |. 

By using the Laplace transform (j20l) and the equation (|22l) , we obtain the non-linear reaction-subdiffusion equation: 



dp 
dt 



D (7) e^o rip(x^u))dujjl^J_^p ^) g- /J ripix,u))du _^ ^ 



(23) 



where D (7) = a'^ /2tq is the anomalous difFusivity. 
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III. MODEL B 



Model B deals with the case when the transport process described by CTRW depends on the chemical reaction. 
This model was considered by Vlad and Ross . They introduced the notion of the age of the particle as the transition 
time between two successive jumps. The particles have zero age when they just arrive at some point x from which 
they will jump later. The main assumption for Model B is that the newborn particles produced with the rate r+ (p) p 
have zero age. In other words, when a new particle is born, it is given a new waiting time for a jump (zero age). 
The density j {x,t) for Model B can be interpreted as a zero- age density of particles arriving at the point x exactly 
at time t. The particles arrive at the point x because of the jumps in space and a birth process with the rate r+ (p) . 
Since the non-linear reaction rate r (p) is the difference between the birth rate r"*" (p) and the death rate (p) , the 
balance equations for the densities j (x, t) and p {x, t) can be written as 

p{x,t)= pa{x)e-^»'' (t) + / j{x,T)e-^> ("(^■""''"l' (t ~ r) dr (24) 

and 







j{x,t) = r+{p)p + J po{x-z)e-^o^ (''("■"»'*"u;(z)0 (t) 

+ /* / jix-z,T)e- (''("■"»'^"u;(z)0 {t - r) dzdr. (25) 



This system of equations has been derived in [3, [T^l.The authors employed the Markov model with age-dependent 
density ^ (x, t, r) such that p (x, t) = ^ (x, t, r) dr. Of course, ([M)) and (^5)) can be formulated directly as the 
balance equations without introduction of ^(x,i, t). Eq. (|24|) is the conservation law for the particles at point x 
at time t. Eq. (|25p describes the situation when the particles with zero age at point x are either produced with a 
rate r+ (p) or arrive at point x from other positions. By using the method developed for Model A, one can derive the 
non-linear Master equation for the density p (x, t) 







p{x-Z,T)e- r-{p{x-z.u))du^ _ ^ g- r- {p{x^u))du 



dr 



+r+ (p) p - r- (p) p. (26) 



Note that Vlad and Ross stated that the balance equations (|24|) and (|25|) could not be reduced to a non-linear Master 

equation as (j26p due to the non-linear term exp^— J^r~ {p (x ~ z,u)) duj [2]. However, Yadav and Horsthemke 

managed to overcome this difficulty by using the large-spatial scale and long-time limits [iB] • They used the standard 
asymptotics for the Fourier transform of jump density w (k) = 1 — a^k^ + o{k'^) and neglected the initial conditions 
in the long-time limit. If we expand p{x ~ z,t) in (j26p for small z, then we obtain the equation derived in [T6| 

§ = Y ^ ^* if (i - r) p (x, r) e- '--(P(-.«))'^-dr + r+ {p) p - (p) p. (27) 

Thus, the Master equation (|26p can be considered as the generalization of (P7)) derived in [16] for the arbitrary jump 
distribution w (z). 

The essence of the Model B and the main difference with Model A is that it describes the situation when the 
newborn particles are given the new waiting time for a jump. One can think of a situation in which the trapping 
mechanism has been induced by chemical binding of newborn molecules. In this case we have to take into account 
the aging effects. Of course, we should make a clear distinction between the age of a jump event and the age of a 
particle from t = 0. The latter effect is not considered in this paper. 



IV. MODEL C 



This model corresponds to the phenonienological generalization of the CTRW model for the case when the chemical 
reaction is taken into account. We can incorporate the local growth rate of diffusing particles by adding a term 
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/o ^'^ (P i^^'^ ~ ^)) P (2^1 i — t) ^' (r) dr to the balance equation for the density of particles. We write 



p{x,t) — p {x,0) 'i' (t) + / / p (x — z,t — t) w{z)(j> (t) dzdr 

Jq Jr 

+ [ r+{p{x,t-T))p{x,t-T)'^{T)dT. (28) 



This equation has been used in [^,0] as a starting point. Henry, Langlands and Wearne modified it in [Tsj and pointed 
out that the Model C is not justified at the mesoscopic level and what is more its interpretation is not clear. The 
purpose of this section is to show that Model C has a very natural physical interpretation. 

Here we assume that the reaction is a pure birth process: r"*" (p) , and j {x,t) is a zero-age density of particles 
arriving at the point x exactly at time t. The balance equations are 

p (x, t) = po {x) ^{t)+ j 3 {x, i - r) * (r) dr (29) 
Jo 

and 

j{x,t) = r+{p)p+ / po (x - z) w{z)(f> (t) dz 
Jr 

t 

j (x — z,t — t) w{z)(j) (t) dzdr. (30) 

These equations do not involve non-linear terms inside the integrals. Therefore, the standard technique of Fourier- 
Laplace transforms can be employed to reduce the two balance equations and (|30p to a single equation for p {x, t). 
It turns out that this equation can be written as a phenomenological balance equation (^5)1 . So Eq. (^5]) corresponds 
to a mesoscopic situation when the newborn particles are given a new waiting time for a jump. Note that it has been 
found [III that we cannot use the balance equation with negative reaction term, since it leads to the negative 
density for the sub-diffusive transport. Of course. Model C is just a particular case of Model B when the death rate 
r-(p) = 0. 

From (|28p one can obtain the Master equation for the density p (x, t) in the foUowing form 
dp ' f 

— = / K (t — t) I p(x — z,t)w {z)dz ~ p(x,t) dr + (p) p 
ot Jo Ur j 

in which the transport term does not directly depend on the chemical reaction as in Models A and B. 



V. SPEED OF TRAVELING WAVES 



In this section we address the problem of wavefront propagation for Model A. We assume that the reaction rate 
r (p) is of the Kolmogorov-Petrovskii-Piskunov (KPP) type [lol. |36|: 

max r (p) = r (0) > 0, r(l) = 0. (31) 

0<p<l 

Note that the standard logistic growth corresponds to r (p) = 1 — p [s^]- We start with the non-linear Master equation 
for the density p(x,t) 

§ = Y ^ ^* ^ (t - r) p (x, r) efr r(p(^:u))dudr + r (p) p (32) 

with an initial condition in the form of a step function 

p(x,0)=^(x), (33) 

where 9 {x) = 1 for x < and 9{x) = for x > 0. This condition describes the initial segregation of an unstable state 
(p = 0) for a; > and a stable state (p = 1) for a; < 0. 
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The purpose is to find the travehng wave solution p{x,t) = f {x — vt) of the initial value problem (l32l) and p3|) . 
Here v is the speed at which the wave profile / invades the unstable state with p — 0. When the memory kernel 
K {t — t) = TQ^5{t — r), then Eq. ([32]) becomes the KPP-equation (Fisher equation) 

|^Z,g..„„ (3., 

with the diffusion coefficient D = (T^/2ro. This case corresponds to the exponential waiting time pdf (i) = A exp (—At) 
with A = Tj^^. It is well-known that the minimal propagation speed v for (|34p is 2\J Dr (0). 

To find the propagation rate for the non-Markovian initial value problem (15^ and (155)) , we use the Hamilton- Jacobi 
approach [H, [sl]. The starting point is to apply hyperbolic scaling x xje, t — !> tje with e ^ which corresponds 
to the long-time large-distance behavior of the traveling wave. When the scaling parameter e — >■ 0, the rescaled 
density p'^{x,t) = pixje^tje) can take only two values and 1 everywhere except in the narrow front region where 
the transport and reaction terms are balanced. In another words, the wave profile / ((x — vt) / e) tends to a unit step 
function 9{x — vt). The aim is to find the location of the front and rate at which it moves. 

We introduce the action functional as 



p'^{x, t) — exp 



G^{x,t) 



(35) 



It follows from ([55)1 that if the function G{x,t) = lim^^o G'^ {x,t) is positive, the rescaled density p'^{x,t) as 
e — >■ 0. The boundary of the set G{x, t) > can be regarded as the reaction front [H, [s^. Then, the front position 
x{t) can be determined by the equation G{x{t),t) = 0. Since we are interested in the leading edge of the traveling 
wave {p^ ~ 0), we write the equation for the rescaled density in the linear form 

Substituting (|35l) into (|36|) and taking the limit e — !> 0, we obtain 



This is the Hamilton- Jacobi equation for the action functional G{x, t) (see [1,[H,[2^). If we introduce the Hamiltonian 
H — and the generalized momentum P = ^ , then equation ([57]) can be written as 

2 2 

H + k{H-r{Q))'^-^-r{Q) = Q, (38) 

where K (s) is the Laplace transform of the memory kernel K (t) . This equation is different from the analogous one 
in 1,111. 

The propagation rate v can be found from [23| 



dH _H 
dp p 

Solving quadratic equation (|38l) for p, we obtain from ([M)) that 



(39) 



where H* is the solution of 



V = i^V J^i^:^, (40) 



d ( H-riO) \ _ 2{H-riO)) ^^^^ 



dH \K{H-r) (0) / HK{H - r(0)) 
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A. Front propagation rate for subdiffusive transport 

Let us consider the reaction-subdiffusion case for which the Laplace transform of the waiting time pdf (j) (t) is 

Hs)= ^^} .y , 0<7<1 

1 + (stq) 

and the Laplace transform of the memory kernel is given by ([20l) . Then Eq. (|4T|) has the solution H* — 2r(0)/ (2 — 7) . 
Substitution of this solution into (PU)) gives the propagation rate v corresponding to the Master equation (15^ with 
the initial condition ([55]) . We introduce the following notations for the propagation speeds: va for Model A; vb and 
vc for Model B and Model C respectively. The minimal propagation speed va for Model A with (|32p and (|33| is 



2r(0)2-7a2 



For simplicity, we consider the case when the death rate (p) obeys r~ (0) — 0. For the Schlogl first model (jS)), it 
means that ^3 = (see (g]). Then for Model B and Model C we have 0, 0] 



/r(0)2-7a2 (3-7)^" 



The case ^ was considered in psl j. 

For 7 = 1, we have the classical result u — 2^y Dr(0) with = a"^ /2tq that corresponds to the KPP-equation 
(Fisher equation) (|34|) . For < 7 < 1, the propagation speed vb is greater than u^. This is because the newborn 
particles in Model B and Model C are given new waiting times for the jump event. As a result the overall transport 
process and the propagation rate for Model A are slower than those of Model B and Model C. It is also interesting 
to compare our results with those obtained in poj . The authors considered the irreversible autocatalytic reaction 
A-\- X ^ 2X with subdiffusion. They found that the minimal propagation is zero. This finding seems to contradict 
with our result of finite speed propagation (|42|) . In fact, Model A and the subdiffusion-reaction model studied in 
[29| are different. In our model we keep the concentration of one of the component constant, while Froemberg et 
al considered two-component system of equations for which both reactants A and X vary in space and time. For 
example, in our paper the density pA of the catalyst A in ([3]) is assumed to be constant. 



VI. CONCLUSIONS 



In this paper we have given a mesoscopic description of a reaction-transport system of particles performing a 
continuous time random walk (CTRW) with non-exponential waiting time distributions. Our main objective has been 
to implement a non-linear kinetic term into non-Markovian transport equations. We have considered three different 
CTRW models with reactions which have been discussed in the literature. We have derived new non-linear Master 
equations (fT4| and (j26p for the mesoscopic density of reacting particles corresponding to CTRW with arbitrary jump 
and waiting time distributions. We have applied the theory to the problem of front propagation in reaction-transport 
systems with KPP-kinetics and non-Markovian diffusion. We have found an explicit expression for the speed of the 
propagating front in the case of subdiffusion transport. 
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